home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / curvatur.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-14  |  31.1 KB  |  772 lines
  1. /******************************************************************************
  2. * Curvature.c - curvature computation of curves and surfaces.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8.  
  9. #define MAX_POS_REF_ITERATION 20
  10.  
  11. /******************************************************************************
  12. * DESCRIPTION:                                                               M
  13. * Computes a scalar curve representing the curvature of a planar curve.      M
  14. *   The given curve is assumed to be planar and only its x and y coordinates M
  15. * are considered.                                 M
  16. *   Then the curvature k is equal to                         M
  17. *      .  ..    .  ..                                 V
  18. *      X  Y  -  Y  X                                 V
  19. * k =  -------------                                 V
  20. *         .2   .2  3/2                                  V
  21. *      (  X  + Y  )                                 V
  22. *   Since we cannot represent k because of the square root, we compute and   M
  23. * represent k^2.                                 M
  24. *                                                                            *
  25. * PARAMETERS:                                                                M
  26. *   Crv:       To compute the square of teh curvature field for.             M
  27. *                                                                            *
  28. * RETURN VALUE:                                                              M
  29. *   CagdCrvStruct *:   The square of the curvature field of Crv.             M
  30. *                                                                            *
  31. * KEYWORDS:                                                                  M
  32. *   SymbCrv2DCurvatureSqr, curvature                                         M
  33. *****************************************************************************/
  34. CagdCrvStruct *SymbCrv2DCurvatureSqr(CagdCrvStruct *Crv)
  35. {
  36.     CagdBType
  37.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  38.     CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
  39.           *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
  40.           *CTmp1, *CTmp2, *CTmp3, *Numer, *Denom, *CurvatureSqr,
  41.     *Crv1Deriv = CagdCrvDerive(Crv),
  42.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  43.  
  44.     SymbCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
  45.     SymbCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
  46.     CagdCrvFree(Crv1Deriv);
  47.     CagdCrvFree(Crv2Deriv);
  48.  
  49.     CTmp1 = SymbCrvMult(Crv1X, Crv2Y);
  50.     CTmp2 = SymbCrvMult(Crv2X, Crv1Y);
  51.     CTmp3 = SymbCrvSub(CTmp1, CTmp2);
  52.     CagdCrvFree(CTmp1);
  53.     CagdCrvFree(CTmp2);
  54.     Numer = SymbCrvMult(CTmp3, CTmp3);
  55.     CagdCrvFree(CTmp3);
  56.  
  57.     CTmp1 = SymbCrvMult(Crv1X, Crv1X);
  58.     CTmp2 = SymbCrvMult(Crv1Y, Crv1Y);
  59.     CTmp3 = SymbCrvAdd(CTmp1, CTmp2);
  60.     CagdCrvFree(CTmp1);
  61.     CagdCrvFree(CTmp2);
  62.     CTmp1 = SymbCrvMult(CTmp3, CTmp3);
  63.     Denom = SymbCrvMult(CTmp1, CTmp3);
  64.     CagdCrvFree(CTmp1);
  65.     CagdCrvFree(CTmp3);
  66.     if (IsRational) {
  67.     CTmp1 = SymbCrvMult(Crv1W, Crv1W);
  68.     CTmp2 = SymbCrvMult(CTmp1, CTmp1);
  69.     CTmp3 = SymbCrvMult(CTmp2, Numer);
  70.     CagdCrvFree(CTmp1);
  71.     CagdCrvFree(CTmp2);
  72.     CagdCrvFree(Numer);
  73.     Numer = CTmp3;
  74.  
  75.     CTmp1 = SymbCrvMult(Crv2W, Crv2W);
  76.     CTmp2 = SymbCrvMult(CTmp1, Denom);
  77.     CagdCrvFree(CTmp1);
  78.     CagdCrvFree(Denom);
  79.     Denom = CTmp2;
  80.     }
  81.     if (CAGD_IS_BSPLINE_CRV(Denom)) {
  82.     CTmp1 = SymbMakePosCrvCtlPolyPos(Denom);
  83.     CagdCrvFree(Denom);
  84.     Denom = CTmp1;
  85.     }
  86.  
  87.     CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
  88.     CurvatureSqr = SymbCrvMergeScalar(Denom, Numer, NULL, NULL);
  89.     CagdCrvFree(Denom);
  90.     CagdCrvFree(Numer);
  91.  
  92.     CagdCrvFree(Crv1X);
  93.     CagdCrvFree(Crv1Y);
  94.     CagdCrvFree(Crv2X);
  95.     CagdCrvFree(Crv2Y);
  96.     if (Crv1Z)
  97.     CagdCrvFree(Crv1Z);
  98.     if (Crv2Z)
  99.     CagdCrvFree(Crv2Z);
  100.     if (Crv1W)
  101.     CagdCrvFree(Crv1W);
  102.     if (Crv2W)
  103.     CagdCrvFree(Crv2W);
  104.  
  105.     return CurvatureSqr;
  106. }
  107.  
  108. /******************************************************************************
  109. * DESCRIPTION:                                                               M
  110. * Computes a vector field curve representing the curvature of a curve, in    M
  111. * the binormal direction, that is kB, and square it:                 M
  112. *                                         M
  113. *       .   ..                     .   ..    2                 V
  114. *       C x C                   ( C x C )                     V
  115. * kB =  -----          and    k^2 = --------                     V
  116. *         .  3                       .  6                     V
  117. *       | C |                      | C |                     V
  118. *                                                                            *
  119. * PARAMETERS:                                                                M
  120. *   Crv:        To compute vector field curvatrue for,                       M
  121. *                                                                            *
  122. * RETURN VALUE:                                                              M
  123. *   CagdCrvStruct *:   Computed vector field curvature of Crv                M
  124. *                                                                            *
  125. * KEYWORDS:                                                                  M
  126. *   SymbCrv3DCurvatureSqr, curvature                                         M
  127. *****************************************************************************/
  128. CagdCrvStruct *SymbCrv3DCurvatureSqr(CagdCrvStruct *Crv)
  129. {
  130.     CagdCrvStruct
  131.     *CTmp, *CTmp2, *Numer, *Denom, *CurvatureSqr,
  132.     *Crv1Deriv = CagdCrvDerive(Crv),
  133.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  134.  
  135.     CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
  136.     CagdCrvFree(Crv2Deriv);
  137.     Numer = SymbCrvDotProd(CTmp, CTmp);
  138.     CagdCrvFree(CTmp);
  139.  
  140.     CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
  141.     CagdCrvFree(Crv1Deriv);
  142.     CTmp2 = SymbCrvMult(CTmp, CTmp);
  143.     Denom = SymbCrvMult(CTmp, CTmp2);
  144.     CagdCrvFree(CTmp);
  145.     CagdCrvFree(CTmp2);
  146.  
  147.     if (CAGD_IS_RATIONAL_CRV(Denom) || CAGD_IS_RATIONAL_CRV(Numer)) {
  148.     CTmp = SymbCrvInvert(Denom);
  149.     CurvatureSqr = SymbCrvMult(CTmp, Numer);
  150.     CagdCrvFree(CTmp);
  151.     }
  152.     else {
  153.     CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ;
  154.  
  155.     SymbCrvSplitScalar(Numer, &PCrvW, &PCrvX, &PCrvY, &PCrvZ);
  156.     CagdMakeCrvsCompatible(&Denom, &PCrvX, TRUE, TRUE);
  157.     CagdMakeCrvsCompatible(&Denom, &PCrvY, TRUE, TRUE);
  158.     CagdMakeCrvsCompatible(&Denom, &PCrvZ, TRUE, TRUE);
  159.     CurvatureSqr = SymbCrvMergeScalar(Denom, PCrvX, PCrvY, PCrvZ);
  160.     CagdCrvFree(PCrvX);
  161.     CagdCrvFree(PCrvY);
  162.     CagdCrvFree(PCrvZ);
  163.     }
  164.  
  165.     CagdCrvFree(Denom);
  166.     CagdCrvFree(Numer);
  167.  
  168.     return CurvatureSqr;
  169. }
  170.  
  171. /******************************************************************************
  172. * DESCRIPTION:                                                               M
  173. * Computes a vector field curve representing the radius (1/curvature) of a   M
  174. * curve, in the normal direction, that is N / k:                 M
  175. *                                         M
  176. *                   .   ..    .         .  6        .   ..    .     .  2         V
  177. *          k N     (C x C ) x C       | C |      ( (C x C ) x C ) | C |         V
  178. * N / k =  ----- = ------------ . --------- = -----------------------         V
  179. *            2          .  4       .   .. 2           .   .. 2             V
  180. *           k         | C |       (C x C )           (C x C )             V
  181. *                                                                            M
  182. *                                                                            *
  183. * PARAMETERS:                                                                M
  184. *   Crv:        To compute the normal field with radius as magnitude.        M
  185. *                                                                            *
  186. * RETURN VALUE:                                                              M
  187. *   CagdCrvStruct *:   Computed normal field with 1 / k as magnitude.        M
  188. *                                                                            *
  189. * KEYWORDS:                                                                  M
  190. *   SymbCrv3DRadiusNormal, curvature                                         M
  191. *****************************************************************************/
  192. CagdCrvStruct *SymbCrv3DRadiusNormal(CagdCrvStruct *Crv)
  193. {
  194.     CagdCrvStruct
  195.     *PCrvW, *PCrvX, *PCrvY, *PCrvZ,
  196.     *CTmp, *CTmp2, *Numer, *Denom, *RadiusNormal,
  197.     *Crv1Deriv = CagdCrvDerive(Crv),
  198.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  199.  
  200.     CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
  201.     CagdCrvFree(Crv2Deriv);
  202.     Denom = SymbCrvDotProd(CTmp, CTmp);
  203.     CTmp2 = SymbCrvCrossProd(CTmp, Crv1Deriv);
  204.     CagdCrvFree(CTmp);
  205.     CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
  206.     Numer = SymbCrvMultScalar(CTmp2, CTmp);
  207.     CagdCrvFree(CTmp2);
  208.     CagdCrvFree(CTmp);
  209.  
  210.     if (CAGD_IS_RATIONAL_CRV(Denom) || CAGD_IS_RATIONAL_CRV(Numer)) {
  211.     CTmp = SymbCrvInvert(Denom);
  212.     RadiusNormal = SymbCrvMult(CTmp, Numer);
  213.     CagdCrvFree(CTmp);
  214.     }
  215.     else {
  216.     SymbCrvSplitScalar(Numer, &PCrvW, &PCrvX, &PCrvY, &PCrvZ);
  217.     CagdMakeCrvsCompatible(&Denom, &PCrvX, TRUE, TRUE);
  218.     CagdMakeCrvsCompatible(&Denom, &PCrvY, TRUE, TRUE);
  219.     CagdMakeCrvsCompatible(&Denom, &PCrvZ, TRUE, TRUE);
  220.     RadiusNormal = SymbCrvMergeScalar(Denom, PCrvX, PCrvY, PCrvZ);
  221.     CagdCrvFree(PCrvX);
  222.     CagdCrvFree(PCrvY);
  223.     CagdCrvFree(PCrvZ);
  224.     }
  225.  
  226.     CagdCrvFree(Denom);
  227.     CagdCrvFree(Numer);
  228.  
  229.     return RadiusNormal;
  230. }
  231.  
  232. /*****************************************************************************
  233. * DESCRIPTION:                                                               M
  234. * Computes a vector field curve representing the curvature of a curve, in    M
  235. * the normal direction, that is kN.                         M
  236. *                .   ..      .       .   ..     .                 V
  237. *                C x C       C     ( C x C  ) x C                      V
  238. * kN = kB x T =  -----  x  ----- = --------------                 V
  239. *                  .  3      .           .  4                     V
  240. *                | C |     | C |       | C |                     V
  241. *                                                                            *
  242. * PARAMETERS:                                                                M
  243. *   Crv:        To compute the normal curvature field.                       M
  244. *                                                                            *
  245. * RETURN VALUE:                                                              M
  246. *   CagdCrvStruct *:   Computed normal curvature field.                      M
  247. *                                                                            *
  248. * KEYWORDS:                                                                  M
  249. *   SymbCrv3DCurvatureNormal, curvature                                      M
  250. *****************************************************************************/
  251. CagdCrvStruct *SymbCrv3DCurvatureNormal(CagdCrvStruct *Crv)
  252. {
  253.     CagdBType
  254.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  255.     CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ,
  256.     *CTmp, *Numer, *Denom, *CurvatureNormal,
  257.     *Crv1Deriv = CagdCrvDerive(Crv),
  258.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  259.  
  260.     CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
  261.     CagdCrvFree(Crv2Deriv);
  262.     Numer = SymbCrvCrossProd(CTmp, Crv1Deriv);
  263.     CagdCrvFree(CTmp);
  264.     SymbCrvSplitScalar(Numer, &CrvW, &CrvX, &CrvY, &CrvZ);
  265.     CagdCrvFree(Numer);
  266.  
  267.     CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
  268.     CagdCrvFree(Crv1Deriv);
  269.     Denom = SymbCrvMult(CTmp, CTmp);
  270.  
  271.     if (IsRational) {
  272.     CagdCrvStruct *DenomCrvW, *DenomCrvX, *DenomCrvY, *DenomCrvZ;
  273.  
  274.     SymbCrvSplitScalar(Denom,
  275.                &DenomCrvW, &DenomCrvX, &DenomCrvY, &DenomCrvZ);
  276.     CagdCrvFree(Denom);
  277.  
  278.     CTmp = SymbCrvMult(CrvW, DenomCrvX);
  279.     CagdCrvFree(CrvW);
  280.     CrvW = CTmp;
  281.  
  282.     CTmp = SymbCrvMult(CrvX, DenomCrvW);
  283.     CagdCrvFree(CrvX);
  284.     CrvX = CTmp;
  285.  
  286.     CTmp = SymbCrvMult(CrvY, DenomCrvW);
  287.     CagdCrvFree(CrvY);
  288.     CrvY = CTmp;
  289.  
  290.     CTmp = SymbCrvMult(CrvZ, DenomCrvW);
  291.     CagdCrvFree(CrvZ);
  292.     CrvZ = CTmp;
  293.  
  294.     CagdCrvFree(DenomCrvW);
  295.     CagdCrvFree(DenomCrvX);
  296.  
  297.     CurvatureNormal = SymbCrvMult(CTmp, Numer);
  298.     CagdCrvFree(CTmp);
  299.     }
  300.     else {
  301.     CagdMakeCrvsCompatible(&Denom, &CrvX, TRUE, TRUE);
  302.     CagdMakeCrvsCompatible(&Denom, &CrvY, TRUE, TRUE);
  303.     CagdMakeCrvsCompatible(&Denom, &CrvZ, TRUE, TRUE);
  304.     CrvW = Denom;
  305.     }
  306.     CurvatureNormal = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);
  307.     CagdCrvFree(CrvX);
  308.     CagdCrvFree(CrvY);
  309.     CagdCrvFree(CrvZ);
  310.     CagdCrvFree(CrvW);
  311.  
  312.     return CurvatureNormal;
  313. }
  314.  
  315. /******************************************************************************
  316. * DESCRIPTION:                                                               M
  317. * Computes a scalar curve representing the curvature sign of a planar curve. M
  318. *   The given curve is assumed to be planar and only its x and y coordinates M
  319. * are considered.                                 M
  320. *   Then the curvature sign is equal to                         M
  321. *      .  ..    .  ..                                 V
  322. * s =  X  Y  -  Y  X                                 V
  323. *                                                                            *
  324. * PARAMETERS:                                                                M
  325. *   Crv:                                                                     M
  326. *                                                                            *
  327. * RETURN VALUE:                                                              M
  328. *   CagdCrvStruct *:                                                         M
  329. *                                                                            *
  330. * KEYWORDS:                                                                  M
  331. *   SymbCrv2DCurvatureSign, curvature                                        M
  332. *****************************************************************************/
  333. CagdCrvStruct *SymbCrv2DCurvatureSign(CagdCrvStruct *Crv)
  334. {
  335.     CagdBType
  336.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  337.     CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
  338.           *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
  339.           *CTmp1, *CTmp2, *Numer, *Denom, *CurvatureSign,
  340.     *Crv1Deriv = CagdCrvDerive(Crv),
  341.     *Crv2Deriv = CagdCrvDerive(Crv1Deriv);
  342.  
  343.     SymbCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
  344.     SymbCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
  345.     CagdCrvFree(Crv1Deriv);
  346.     CagdCrvFree(Crv2Deriv);
  347.  
  348.     CTmp1 = SymbCrvMult(Crv1X, Crv2Y);
  349.     CTmp2 = SymbCrvMult(Crv2X, Crv1Y);
  350.     Numer = SymbCrvSub(CTmp1, CTmp2);
  351.     CagdCrvFree(CTmp1);
  352.     CagdCrvFree(CTmp2);
  353.  
  354.     if (IsRational) {
  355.     Denom = SymbCrvMult(Crv1W, Crv2W);
  356.  
  357.     CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
  358.     CurvatureSign = SymbCrvMergeScalar(Denom, Numer, NULL, NULL);
  359.     CagdCrvFree(Denom);
  360.     CagdCrvFree(Numer);
  361.     }
  362.     else {
  363.     CurvatureSign = Numer;
  364.     }
  365.  
  366.     CagdCrvFree(Crv1X);
  367.     CagdCrvFree(Crv1Y);
  368.     CagdCrvFree(Crv2X);
  369.     CagdCrvFree(Crv2Y);
  370.     if (Crv1Z)
  371.     CagdCrvFree(Crv1Z);
  372.     if (Crv2Z)
  373.     CagdCrvFree(Crv2Z);
  374.     if (Crv1W)
  375.     CagdCrvFree(Crv1W);
  376.     if (Crv2W)
  377.     CagdCrvFree(Crv2W);
  378.  
  379.     return CurvatureSign;
  380. }
  381.  
  382. /******************************************************************************
  383. * DESCRIPTION:                                                               M
  384. * Given a scalar curve that is positive, refine it until all its control     M
  385. * points has positive coefficients. Always returns a Bspline curve.          M
  386. *                                                                            *
  387. * PARAMETERS:                                                                M
  388. *   OrigCrv:    To refine until all its control points are non negative.     M
  389. *                                                                            *
  390. * RETURN VALUE:                                                              M
  391. *   CagdCrvStruct *:    Refined positive curve with positive control points. M
  392. *                                                                            *
  393. * KEYWORDS:                                                                  M
  394. *   SymbMakePosCrvCtlPolyPos, refinement                                     M
  395. *****************************************************************************/
  396. CagdCrvStruct *SymbMakePosCrvCtlPolyPos(CagdCrvStruct *OrigCrv)
  397. {
  398.     int l;
  399.     CagdCrvStruct
  400.     *RefCrv = NULL;
  401.  
  402.     switch (OrigCrv -> GType) {
  403.     case CAGD_CBEZIER_TYPE:
  404.         RefCrv = CnvrtBezier2BsplineCrv(OrigCrv);
  405.         break;
  406.     case CAGD_CBSPLINE_TYPE:
  407.         RefCrv = CagdCrvCopy(OrigCrv);
  408.         break;
  409.     case CAGD_CPOWER_TYPE:
  410.     default:
  411.         SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
  412.         break;
  413.     }
  414.  
  415.     for (l = 0; l < MAX_POS_REF_ITERATION; l++) {
  416.     int i, j,
  417.         Len = RefCrv -> Length,
  418.         Order = RefCrv -> Order;
  419.     CagdRType
  420.         *KV = RefCrv -> KnotVector,
  421.         *Nodes = BspKnotNodes(KV, Len + Order, Order),
  422.         *Pts = RefCrv -> Points[1];
  423.  
  424.     for (i = j = 0; i < Len; i++) {
  425.         if (ABS(Pts[i]) < SQR(EPSILON))
  426.         Pts[i] = 0.0;
  427.         if (Pts[i] < 0)        /* To refine at negative control points. */
  428.             Nodes[j++] = Nodes[i];
  429.     }
  430.     if (j == 0) {
  431.         IritFree(Nodes);
  432.         break;
  433.     }
  434.     else {
  435.         CagdCrvStruct
  436.             *NewCrv = CagdCrvRefineAtParams(RefCrv, FALSE, Nodes, j);
  437.  
  438.         CagdCrvFree(RefCrv);
  439.         RefCrv = NewCrv;
  440.         IritFree(Nodes);
  441.     }
  442.     }
  443.  
  444.     return RefCrv;
  445. }
  446.  
  447. /*****************************************************************************
  448. * DESCRIPTION:                                                               M
  449. * Given a planar curve, finds all its inflection points by finding the zero  M
  450. * set of the sign of the curvature function of the curve.             M
  451. *                                                                            *
  452. * PARAMETERS:                                                                M
  453. *   Crv:        To find all its inflection points.                           M
  454. *   Epsilon:    Accuracy control.                                            M
  455. *                                                                            *
  456. * RETURN VALUE:                                                              M
  457. *   CagdPtStruct *:    A list of parameter values on Crv that are inflection M
  458. *                      points.                                               M
  459. *                                                                            *
  460. * KEYWORDS:                                                                  M
  461. *   SymbCrv2DInflectionPts, curvature, inflection points                     M
  462. *****************************************************************************/
  463. CagdPtStruct *SymbCrv2DInflectionPts(CagdCrvStruct *Crv, CagdRType Epsilon)
  464. {
  465.     CagdCrvStruct
  466.     *CrvtrSign = SymbCrv2DCurvatureSign(Crv);
  467.     CagdPtStruct
  468.     *InflectionPts = SymbCrvZeroSet(CrvtrSign, 1, Epsilon);
  469.  
  470.     CagdCrvFree(CrvtrSign);
  471.  
  472.     return InflectionPts;    
  473. }
  474.  
  475. /*****************************************************************************
  476. * DESCRIPTION:                                                               M
  477. * Given a planar curve, finds all its extreme curvatrue points by finding    M
  478. * the set of extreme locations on  the curvature function of Crv.         M
  479. *                                                                            *
  480. * PARAMETERS:                                                                M
  481. *   Crv:       To find all int extrem curvature locations.                   M
  482. *   Epsilon:   Accuracy control.                                             M
  483. *                                                                            *
  484. * RETURN VALUE:                                                              M
  485. *   CagdPtStruct *:   A list of parameter values on Crv that have extrem     M
  486. *                     curvature values.                                      M
  487. *                                                                            *
  488. * KEYWORDS:                                                                  M
  489. *   SymbCrv2DExtremCrvtrPts, curvature                                       M
  490. *****************************************************************************/
  491. CagdPtStruct *SymbCrv2DExtremCrvtrPts(CagdCrvStruct *Crv, CagdRType Epsilon)
  492. {
  493.     if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
  494.     CagdCrvStruct
  495.         *CrvtrSqrCrv = SymbCrv2DCurvatureSqr(Crv);
  496.     CagdPtStruct
  497.         *ExtremCrvtrPts = SymbCrvExtremSet(CrvtrSqrCrv, 1, Epsilon);
  498.  
  499.     CagdCrvFree(CrvtrSqrCrv);
  500.  
  501.     return ExtremCrvtrPts;    
  502.     }
  503.     else if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 3) {
  504.     CagdCrvStruct
  505.         *CrvtrNormalCrv = SymbCrv3DCurvatureNormal(Crv),
  506.         *CrvtrMag = SymbCrvDotProd(CrvtrNormalCrv, CrvtrNormalCrv);
  507.     CagdPtStruct
  508.         *ExtremCrvtrPts = SymbCrvExtremSet(CrvtrMag, 1, Epsilon);
  509.  
  510.     CagdCrvFree(CrvtrNormalCrv);
  511.     CagdCrvFree(CrvtrMag);
  512.  
  513.     return ExtremCrvtrPts;    
  514.     }
  515.     else {
  516.     SYMB_FATAL_ERROR(SYMB_ERR_ONLY_2D_OR_3D);
  517.     return NULL;
  518.     }
  519. }
  520.  
  521. /*****************************************************************************
  522. * DESCRIPTION:                                                               M
  523. * Computes coefficients of the first fundamental form of given surface Srf.  M
  524. *                                                                            *
  525. * PARAMETERS:                                                                M
  526. *   Srf:         Do compute the coefficients of the FFF for.                 M
  527. *   DuSrf:       First derivative of Srf with respect to U goes to here.     M
  528. *   DvSrf:       First derivative of Srf with respect to V goes to here.     M
  529. *   FffG11:      FFF G11 scalar field.                                       M
  530. *   FffG12:      FFF G12 scalar field.                                       M
  531. *   FffG22:      FFF G22 scalar field.                                       M
  532. *                                                                            *
  533. * RETURN VALUE:                                                              M
  534. *   void                                                                     M
  535. *                                                                            *
  536. * KEYWORDS:                                                                  M
  537. *   SymbSrfFff, first fundamental form                                       M
  538. *****************************************************************************/
  539. void SymbSrfFff(CagdSrfStruct *Srf,
  540.         CagdSrfStruct **DuSrf,
  541.         CagdSrfStruct **DvSrf,
  542.         CagdSrfStruct **FffG11,
  543.         CagdSrfStruct **FffG12,
  544.         CagdSrfStruct **FffG22)
  545. {
  546.     *DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR);
  547.     *DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR);
  548.  
  549.     *FffG11 = SymbSrfDotProd(*DuSrf, *DuSrf);
  550.     *FffG12 = SymbSrfDotProd(*DuSrf, *DvSrf);
  551.     *FffG22 = SymbSrfDotProd(*DvSrf, *DvSrf);
  552. }
  553.  
  554. /*****************************************************************************
  555. * DESCRIPTION:                                                               M
  556. * Computes coefficients of the first fundamental form of given surface Srf.  M
  557. *   These coefficients are using non normalized normal that is also returned.M
  558. *                                                                            *
  559. * PARAMETERS:                                                                M
  560. *   DuSrf:       First derivative of Srf with respect to U.             M
  561. *   DvSrf:       First derivative of Srf with respect to V.             M
  562. *   SffL11:      SFF L11 scalar field returned herein.                       M
  563. *   SffL12:      SFF L12 scalar field returned herein.                       M
  564. *   SffL22:      SFF L22 scalar field returned herein.                       M
  565. *   SNormal:     Unnormalized normal vector field returned herein.           M
  566. *                                                                            *
  567. * RETURN VALUE:                                                              M
  568. *   void                                                                     M
  569. *                                                                            *
  570. * KEYWORDS:                                                                  M
  571. *   SymbSrfFff, second fundamental form                                      M
  572. *****************************************************************************/
  573. void SymbSrfSff(CagdSrfStruct *DuSrf,
  574.         CagdSrfStruct *DvSrf,
  575.         CagdSrfStruct **SffL11,
  576.         CagdSrfStruct **SffL12,
  577.         CagdSrfStruct **SffL22,
  578.         CagdSrfStruct **SNormal)
  579. {
  580.     CagdSrfStruct
  581.     *DuuSrf = CagdSrfDerive(DuSrf, CAGD_CONST_U_DIR),
  582.     *DuvSrf = CagdSrfDerive(DuSrf, CAGD_CONST_V_DIR),
  583.     *DvvSrf = CagdSrfDerive(DvSrf, CAGD_CONST_V_DIR);
  584.  
  585.     *SNormal = SymbSrfCrossProd(DvSrf, DuSrf);
  586.     *SffL11 = SymbSrfDotProd(DuuSrf, *SNormal);
  587.     *SffL12 = SymbSrfDotProd(DuvSrf, *SNormal);
  588.     *SffL22 = SymbSrfDotProd(DvvSrf, *SNormal);
  589.  
  590.     CagdSrfFree(DuuSrf);
  591.     CagdSrfFree(DuvSrf);
  592.     CagdSrfFree(DvvSrf);
  593. }
  594.  
  595. /*****************************************************************************
  596. * DESCRIPTION:                                                               M
  597. * Computes the expression of Srf11 * Srf22 - Srf12 * Srf21, which is a         M
  598. * determinant of a 2 by 2 matrix.                         M
  599. *                                                                            *
  600. * PARAMETERS:                                                                M
  601. *   Srf11, Srf12, Srf21, Srf22:  The four factors of the determinant.        M
  602. *                                                                            *
  603. * RETURN VALUE:                                                              M
  604. *   CagdSrfStruct *: A scalar field representing the determinant computation.M
  605. *                                                                            *
  606. * KEYWORDS:                                                                  M
  607. *   SymbSrf2DDeterminant, determinant                                        M
  608. *****************************************************************************/
  609. CagdSrfStruct *SymbSrf2DDeterminant(CagdSrfStruct *Srf11,
  610.                     CagdSrfStruct *Srf12,
  611.                     CagdSrfStruct *Srf21,
  612.                     CagdSrfStruct *Srf22)
  613. {
  614.     CagdSrfStruct
  615.     *Prod1 = SymbSrfMult(Srf11, Srf22),
  616.     *Prod2 = SymbSrfMult(Srf21, Srf12),
  617.     *Add12 = SymbSrfSub(Prod1, Prod2);
  618.  
  619.     CagdSrfFree(Prod1);
  620.     CagdSrfFree(Prod2);
  621.     return Add12;
  622. }
  623.  
  624. /*****************************************************************************
  625. * DESCRIPTION:                                                               M
  626. * Computes curvature upper bound as Xi = k1^2 + k2^2, where k1 and k2 are    M
  627. * the principal curvatures.                              M
  628. * Gij are the coefficients of the first fundamental form and Lij are of the  M
  629. * second, using non unit normal n,                         M
  630. *                                         M
  631. *      ( G11 L22 + G22 L11 - 2 G12 L12 )^2 - 2 |G| |L|                 V
  632. * Xi = -----------------------------------------------                 V
  633. *                 |G|^2 ||n||^2                     V
  634. *                                         M
  635. * See: "Second Order Surface Analysis Using Hybrid of Symbolic and Numeric   M
  636. * Operators", By Gershon Elber and Elaine Cohen, Transaction on graphics,    M
  637. * Vol. 12, No. 2, pp 160-178, April 1993.                     M
  638. *                                                                            *
  639. * PARAMETERS:                                                                M
  640. *   Srf:       Surface to compute curvature bound for.                       M
  641. *                                                                            *
  642. * RETURN VALUE:                                                              M
  643. *   CagdSrfStruct *:   A scalar field representing the curvature bound.      M
  644. *                                                                            *
  645. * KEYWORDS:                                                                  M
  646. *   SymbSrfCurvatureUpperBound, curvature                                    M
  647. *****************************************************************************/
  648. CagdSrfStruct *SymbSrfCurvatureUpperBound(CagdSrfStruct *Srf)
  649. {
  650.     CagdSrfStruct *DuSrf, *DvSrf, *SNormal, *STmp1, *STmp2, *STmp3, *STmp4,
  651.     *Numer, *FffDeterminant, *SffDeterminant, *CurvatureBound,
  652.     *Denom, *FffG11, *FffG12, *FffG22, *SffL11, *SffL12, *SffL22;
  653.  
  654.     SymbSrfFff(Srf, &DuSrf, &DvSrf, &FffG11, &FffG12, &FffG22);
  655.     SymbSrfSff(DuSrf, DvSrf, &SffL11, &SffL12, &SffL22, &SNormal);
  656.     CagdSrfFree(DuSrf);
  657.     CagdSrfFree(DvSrf);
  658.  
  659.     STmp1 = SymbSrfMult(FffG11, SffL22);
  660.     STmp2 = SymbSrfMult(FffG22, SffL11);
  661.     STmp3 = SymbSrfMult(FffG12, SffL12);
  662.     STmp4 = SymbSrfScalarScale(STmp3, 2.0);
  663.     CagdSrfFree(STmp3);
  664.     STmp3 = SymbSrfAdd(STmp1, STmp2);
  665.     CagdSrfFree(STmp1);
  666.     CagdSrfFree(STmp2);
  667.     STmp1 = SymbSrfSub(STmp3, STmp4);
  668.     CagdSrfFree(STmp3);
  669.     CagdSrfFree(STmp4);
  670.     STmp2 = SymbSrfMult(STmp1, STmp1);
  671.     CagdSrfFree(STmp1);
  672.  
  673.     FffDeterminant = SymbSrf2DDeterminant(FffG11, FffG12, FffG12, FffG22);
  674.     SffDeterminant = SymbSrf2DDeterminant(SffL11, SffL12, SffL12, SffL22);
  675.     CagdSrfFree(FffG11);
  676.     CagdSrfFree(FffG12);
  677.     CagdSrfFree(FffG22);
  678.     CagdSrfFree(SffL11);
  679.     CagdSrfFree(SffL12);
  680.     CagdSrfFree(SffL22);
  681.  
  682.     STmp1 = SymbSrfMult(FffDeterminant, SffDeterminant);
  683.     STmp3 = SymbSrfScalarScale(STmp1, 2.0);
  684.     CagdSrfFree(STmp1);
  685.     Numer = SymbSrfSub(STmp2, STmp3);
  686.     CagdSrfFree(STmp2);
  687.     CagdSrfFree(STmp3);
  688.  
  689.     STmp1 = SymbSrfDotProd(SNormal, SNormal);
  690.     CagdSrfFree(SNormal);
  691.  
  692.     STmp2 = SymbSrfMult(FffDeterminant, FffDeterminant);
  693.     Denom = SymbSrfMult(STmp1, STmp2);
  694.     CagdSrfFree(STmp1);
  695.     CagdSrfFree(STmp2);
  696.  
  697.     CagdMakeSrfsCompatible(&Denom, &Numer, TRUE, TRUE, TRUE, TRUE);
  698.     CurvatureBound = SymbSrfMergeScalar(Denom, Numer, NULL, NULL);
  699.     CagdSrfFree(Denom);
  700.     CagdSrfFree(Numer);
  701.  
  702.     return CurvatureBound;
  703. }
  704.  
  705. /*****************************************************************************
  706. * DESCRIPTION:                                                               M
  707. * Computes normal curvature bound in given isoparametric direction.          M
  708. *   This turns out to be (L11 . n) / G11 for u and (L22 . n) / G22 for v.    M
  709. *   Herein the square of these equations is computed symbolically and        M
  710. * returned.                                     M
  711. *                                                                            *
  712. * PARAMETERS:                                                                M
  713. *   Srf:      To compute normal curvature in an isoparametric direction Dir. M
  714. *   Dir:      Direction to compute normal curvature. Either U or V.          M
  715. *                                                                            *
  716. * RETURN VALUE:                                                              M
  717. *   CagdSrfStruct *:  A scalar field representing the normal curvature       M
  718. *                     square of Srf in dirction Dir.                         M
  719. *                                                                            *
  720. * KEYWORDS:                                                                  M
  721. *   SymbSrfIsoDirNormalCurvatureBound, curvature                             M
  722. *****************************************************************************/
  723. CagdSrfStruct *SymbSrfIsoDirNormalCurvatureBound(CagdSrfStruct *Srf,
  724.                          CagdSrfDirType Dir)
  725. {
  726.     CagdSrfStruct *DuSrf, *DvSrf, *SNormal, *STmp1, *STmp2, *STmp3, *STmp4,
  727.     *SNormalSize, *FffG11, *FffG12, *FffG22, *SffL11, *SffL12, *SffL22,
  728.     *CurvatureBound;
  729.  
  730.     SymbSrfFff(Srf, &DuSrf, &DvSrf, &FffG11, &FffG12, &FffG22);
  731.     SymbSrfSff(DuSrf, DvSrf, &SffL11, &SffL12, &SffL22, &SNormal);
  732.     CagdSrfFree(DuSrf);
  733.     CagdSrfFree(DvSrf);
  734.  
  735.     SNormalSize = SymbSrfDotProd(SNormal, SNormal);
  736.  
  737.     switch (Dir) {
  738.     case CAGD_CONST_U_DIR:
  739.         STmp2 = SymbSrfMult(SffL11, SffL11);
  740.         STmp3 = SymbSrfMult(FffG11, FffG11);
  741.         STmp4 = SymbSrfMult(SNormalSize, STmp3);
  742.         CagdSrfFree(STmp3);
  743.         STmp1 = SymbSrfInvert(STmp4);
  744.         CagdSrfFree(STmp4);
  745.         break;
  746.     case CAGD_CONST_V_DIR:
  747.         STmp2 = SymbSrfMult(SffL22, SffL22);
  748.         STmp3 = SymbSrfMult(FffG22, FffG22);
  749.         STmp4 = SymbSrfMult(SNormalSize, STmp3);
  750.         CagdSrfFree(STmp3);
  751.         STmp1 = SymbSrfInvert(STmp4);
  752.         CagdSrfFree(STmp4);
  753.         break;
  754.     default:
  755.         SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
  756.         STmp1 = STmp2 = NULL;
  757.     }
  758.  
  759.     CurvatureBound = SymbSrfMult(STmp1, STmp2);
  760.     CagdSrfFree(STmp1);
  761.     CagdSrfFree(STmp2);
  762.     CagdSrfFree(SNormalSize);
  763.     CagdSrfFree(FffG11);
  764.     CagdSrfFree(FffG12);
  765.     CagdSrfFree(FffG22);
  766.     CagdSrfFree(SffL11);
  767.     CagdSrfFree(SffL12);
  768.     CagdSrfFree(SffL22);
  769.  
  770.     return CurvatureBound;
  771. }
  772.